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ABSTRACT 


The  following  report  describes  (a)  an  outline  of  the  Kalman-based 
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cf  a  mathematical  model  to  implement  this  procedure  in  radar  polar 
coordinates,  and  (c)  methods  for  checking,  via  simulation,  the  degree 
of  degradation  introduced  by  model  simplifications  and  other  changes. 
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I.  RECURSIVE  TRACKING 


1.  Introduction 

The  following  report  describes  the  derivation  and  implementation  of 
methods  and  models  for  simulating  the  reentry  tracking  of  an  object  in 
radar -centered  polar  coordinates  and  the  statistical  comparison  of  various 
tracking  algorithms.  This  work  was  performed  for  Lincoln  Laboratory, 

M*l.  T.,  under  the  direction  of  Dr.  Richard  Wishner,  Assistant  Group 
Leader,  Group  42,  in  support  of  a  program  for  advancing  the  techniques 
of  reentry  tracking  from  radar  data. 

The  target  tracking  procedure  considered  here  is  a  relatively  straight- 

forward  application  of  Kalman's  recursive  estimation  technique'  ' .  This 

(21 

application  was  reported  upon  previously  by  Gruber'  forming  the  basis 

for  subsequent  computer  simulation  models.  The  early  models  were  based 

on  a  radar- centered  rectangular  coordinate  ‘stem  and  were  prograznmed 

(31 

and  investigated  by  Bertolini' 

A  radar-centered  polar  coordinate  system  has  the  relative  advantage  of 
leading  to  a  simple  linear  model  for  observation  geometry  at  the  expense  of 
a  more  complicated  dynamical  model.  It  has  been  suspected  that  the  linear 
observation  model  represents  a  definite  advantage  in  tracking  stability  and 
accuracy,  especially  just  after  target  acquisition.  The  present  work  is  a 
contribution  to  the  polar  coordinate  studies. 

Dynamical  equations  for  a  trajectory  in  radar  polar  coordinates  have 
been  given  by  Catalano^^\  We  rcderive  these  equations  here  in  a  matrix 
format  and  continue  with  a  development  of  the  partial  derivative  matrices 
for  the  spherical  earth  case,  which  are  necessary  for  the  tracking 
algorithms. 

In  stimmary,  the  report  contains  (a)  an  outline  of  the  Kalman-based 
tracking  procedure  for  reentry  trajectories,  (b)  a  detailed  derivation  of 
a  mathematical  model  to  implement  this  procedure  in  radar  polar  coor¬ 
dinates,  (c)  methods  for  checking,  via  sim\ilation,  the  degree  of  degradation 
introduced  by  model  simplifications  and  other  changes. 
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An  Appendix  summarizes  the  rather  extensive  notation  conventions 
and  symbology  used  in  this  report. 

2.  Tracking  Models 

Tile  tracking  procedure,  as  we  noted  above,  is  based  on  Kalman's 
method  of  recursive  estimation.  This  estimation  technique  is  ideally 
suited  to  the  present  problem  because 

a.  It  is  a  general  method.  Many  different  variations  of  reentry 
dynamics  and  observational  procedures  can  be  handled  within 
a  common  framework.  Thus  design  experimentation  is 
facilitated. 

b.  It  is  real-time  estimation  method.  The  utilization  of  tracking 
data  in  forming  the  track  estimate  is  sequential  in  time.  Since 
we  are  interested  in  applying  the  results  of  our  studies  to  real¬ 
time  systems,  we  are  obliged  to  confine  ourselves  to  such 
methods. 

c.  It  provides  optimal  estimates  where  the  model  equations  are 
linear,  and  approximate,  optimal  estimates  in  the  nonlinear 
case. 

In  view  of  these  properties  of  recursive  estimation  it  is  not  surprising 
that  most  of  the  tracking  methods  formulated  in  the  past  for  real-time 
applications  turn  out  to  be  special  cases  of  the  general  Kalman  formulation. 
While  it  is  possible  to  conceive  of  tracking  schemes  which  do  not  fit  such 
a  formulation,  the  investigation  of  each  such  scheme  woizld  be  a  new  and 
unique  undertaking.  We  believe  that  at  present  the  best  service  is  per¬ 
formed  by  investigating  the  Kalman  class  of  tracking  models  and  will 
confine  our  discussion  accordingly. 

In  order  to  implement  a  particular  model,  we  must  postulate  or  derive 
two  sets  of  equations.  The  first  set,  a  system  of  first  order  differential 
equations,  describes  the  dynamics  of  the  reentry  body  and  its  motion  along 
the  trajectory.  The  second  set  describes  the  geometry  of  each  measure¬ 
ment.  Since  it  is  clear  that  the  same  physical  trajectory  can  be  treated  in 
different  coordinate  frameworks  and  that  varying  levels  of  detail  can  be 
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incorporated  in  the  dynamics,  a  class  of  models  arises,  each  treating 
the  same  problem. 

The  dynamical  equations  may  be  written  in  vector  form 
(1.1)  ^  =  f(3^  t) 

where  x  components  are  the  dynamical  or  so-called  state  variables  and 
t  is  time.  Similarly,  the  observation  geometry  may  be  written  in  the 
form 

(1.  2)  h  (3^  t)  +  w  (t) 

where  z  is  the  measurement  vector  at  time  t  and  w  is  an  additive  noise 
vector.  This  noise  vector  is  required  to  be  statistically  uncorrelated  for 
difiJerent  t  and  has  the  assigned  covariance  matrix 

(1.  3)  W(t)  =  cov  (w(t) )  . 

Without  loss  of  generality  w(t)  can  be  assumed  to  have  zero  mean  for 
each  t. 

The  dynaonical  equations  may  be  generalized  to  account  for  random 
disturbances  by  including  random  driving  terms  on  the  right  side. 

(1. 4)  ^ 

In  this  more  general  forix),  ^(t)  is  assumed  to  be  a  zero  mean  white  noise 
process  with  given  covariance  matrix,  Wj^(t).  In  the  following  discussion, 
the  simpler  dynaznical  equations  (1. 1)  will  be  used.  Then  the  modifications 
necessary  to  include  the  effects  of  a  random  driving  terin,  ^^(t),  will  be 
noted. 

Owing  to  the  nonlinear  nature  of  the  above  equations  (f_  and  h  are 
nonlinear  functions  of  x),  it  is  inevitable  that  statistical  calctilations  be 
based  on  a  linearization  of  (1. 1)  and  (1.  2).  Thus,  to  compute  tracking 
errors  and  their  statistics,  we  assume  that  the  reentry  object  exhibits 


-  3  - 


AfsaCOM 


only  small  deviations  from  a  nominal  trajectory  and.  from  (1, 1)  derive 
an  equation  for  the  propagation  of  such  deviations.  Let  us  therefore 
consider  a  specific  trajectory,  ^gCt),  satisfying  (1, 1),  as  a  reference 
trajectory 

(1.  5)  -^  =  f  t) 

The  actual  trajectory  will  differ  from  by  an  amount  62^  assumed  small. 

(1.6)  fix  =  X  -  Xq 

From  (1. 1)  we  have 
dx  , 

"dt  ~  * 

Expanding  terms  on  the  right  in  a  power  series  to  the  first  order  in  fix 
and  using  (1.  5),  we  find  a  matrix  equation  for  the  dynamics  of  fix. 

(1.7)  fix  =  J  (Xq,  t)  fix 


evaluated  on  the  reference  trajectory. 

Similarly,  we  use  as  inputs  to  a  statistical  calculation,  deviations  of 
the  actual  observation  vector,  z,  from  the  ideal  observations,  that 
would  be  taken  on  the  assumed  reference  trajectory  with  no  noise.  That 
is. 


(1.  9)  5  z  =  z  - 

where 
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(1.10)  ^  =  h(^  t)  . 

Again  asstuning  that  6  z  and  6  x  are  small,  making  a  first  order  approxima¬ 
tion  in  (1.  2)  and  using  (1. 6)  and  (1. 10),  we  find  from 

z^  +  6z  =  h(xQ  +  63^  t)  +  w(t) 

that 

(1.11)  bz  =  H(3^,  t)  6x  +  w(^)  • 

Here  H  i£  a  matrix  of  first  partials  of  h  with  components 

(1.12)  = 

evalxiated  on  the  reference  trajectory. 

3.  Tracking  Recursions 

In  outlining  the  method  of  tracking  by  the  recursive  technique  with  the 
above  model,  we  assume  that  data  is  received  at  discrete  times  t^  (n  =  1, 

2,  3,  . . . )  and  introduce  the  following  quantities  and  notation. 

£(n)  -  the  data  received  at  time  n. 

x(n)  -  the  value  of  the  state  vector  at  time  n. 

x(£(n)  -  an  estimate  of  x(n)  given  the  data  up  to  and 

including  £(n). 

x(n'n-l)  -  an  estimate  of  x(n)  given  data  up  to  and  including 
£(n-l). 

2  (njn-l)  -  the  covariance  matrix  of  errors  in  x(n|n-l). 

H(n}  -  partial  derivative  matrix  (1. 12)  at  t^ 

\V(n)  -  measiirement  noise  covariance  matrix  (1.  3}  at  t^. 

t 


In  what  follows. 


denotes  matrix  transposition, 
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We  begin  by  as8\iining  that  we  are  in  possession  of  x(n-l  I  n-1)  and 
S  (n-ljn-1).  Then  the  recursive  method  consists  of  the  following  steps: 

a.  Extrpolate  x  (n-lj  n-1)  forward  to  time  t^  thus  obtaining  x  (n|  n-1). 
The  method  of  extrapolation  is  to  use  x(n-l|n-l)  as  the  initial 
vector  in  the  dynamical  equations  (1. 1)  and  to  numerically  integrate 
these  equations  up  to  time  t^^j^  . 

A 

(1.13)  ^  =  f(x.  t) 

The  solution  at  time  t^^j  is  x  (nl  n-1). 

b.  Extrapolate  S  (n-lj  n-1)  forward  to  time  t^  thus  obtaining 

S  (nln-1).  The  method  here  is  to  adapt  the  linearized  error 
propagation  equation  (1.7).  The  corresponding  propagation 
equation  for  covariances  has  the  form 

(1.14)  ^=J2J^ 

where,  ideally,  J  is  the  Jacobian  (1.  8)  evaluated  on  the  true 
trajectcry.  However,  for  a  real  estimation  probleni,  the  true 
trajectory  is  not  known.  Hence,  we  evaluate  J  instead  on  the 
estimated  trajectory  between  t^  and  t^  determined  by  the 
extrapolation  process  of  a.  Then  E  (n-l|n-l)  is  used  as  an  initial 
n  atrix  value  in  (1. 14)  and  this  matrix  equation  in  numerically 
integrated  forward  to  provide  a  solution  2  (njn-1)  at  time  t^  . 

c.  At  time  t^,  new  data£(n)  are  received;  2(n|n)  andx(n|n)  are 
generated  from  2  (n|  n-1)  and  x (njn-1)  by  the  Kalman  formulas 

(1. 15)  2  (nj  n)  =  2  (nj  n-1)  -  2  (nj  n-1)  H^(n)  [H(n)  2  nl  n-1)  H^n) 

+  W(n)r^  H(n)  2(nln-l) 
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(1. 16)  X  (u]  n)  =  X  (nj  n-1)  +  S  (n|  n-1)  H^n)  [H(ii)  Z  (nj  n-1)  H^n) 

+  W(n)]"^  (z(n)  -  H(n)  x  (nln-1)] 

The  cycle  of  calculations  then  repeats,  returning  to  step  a.  The 
successive  values  of  x{n|n)  thus  generated  are  the  tracking  estixnates. 

Of  course,  an  independent  method  is  required  to  supply  the  initial  starting 
values  x(C  1 0)  and  Z  (0 1  0). 

The  above  procedure  may  be  generalized  slightly  to  account  for  the 
presence  of  random  driving  terms  i:'  the  dynamical  equations  as  in  (1.  4). 

In  most  cases  such  terms  in  the  model  are  added  to  account  for  unknown 
effects  and  to  hedge  the  statistical  estimation  process.  For  such  a  purpose 
it  is  equally  appropriate  to  add  the  random  term  only  at  the  discrete  obser¬ 
vation  times.  Thus  instead  of  analyzing  the  effect  of  a  continuous  excitation 
w^(t),  ve  simply  replace  it  with  a  discrete  one  Wjj.(n)  with  an  appropriate 
covariance  W^{n).  Only  the  above  step  b.  is  modified  to  account  for  this 
complication.  When  Z  (n>l|  n-1)  is  extrapolated  forward,  the  matrix 
W^(n-l)  is  added  to  the  final  extrapolated  value  in  order  to  obtain 
Z  (n|n-l). 
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n.  RADAR  POLAR  COORDINATE  MODEL 
ELLIPSOIDAL  EARTH 


1.  Earth-Station  Geometry 

It  should  be  clear  that  if  the  same  physical  problem  is  described  in 
different  coordinate  systems  the  functions  1^  and  J,  H  of  the  model 
will  differ  in  a  corresponding  way.  Since  the  effects  of  the  linear 
approximations  on  the  accuracy  and  stability  of  subsequent  statistical 
calculations  depend  on  the  precise  notdinear  forms  of  1^  it  is  of 
some  interest  to  explore  various  coordinate  descriptions  for  their 
influence  on  these  numerical  properties.  At  the  same  time,  the  level 
of  complexity  included  in  the  dynamical  model  may  be  varied  to  account 
for  the  various  physical  influences  and  the  properties  of  the  reentry  body. 

In  this  report,  we  shall  describe  in  detail  one  such  coordinate  system 
and  level  of  model  complexity.  The  chosen  coordinate  system  is  earth 
fixed,  polar,  and  centered  at  an  observation  site.  Observations  are 
assumed  tc  be  made  by  a  single  radar  and  consist  of  measurements  of 
range  R,  azimuth  A,  elevation  E,  and  their  rates  (or  a  suitable  subset 
of  these).  By  choosing  these  very  same  variables  as  dynaznical 
variables,  we  insure  that  the  functions  h  are  exactly  linear  and  that 
no  approximation  need  be  made  on  deriving  H  from  ^  On  the  other 
hand,  the  dynamical  laws  (1. 1)  are  complicated  by  this  choice. 

The  dynamical  model  assiunes  a  rotating  earth  of  ellipsoidal  figture, 
gravity  terms  including  the  first  zonal  harmonic  (J^)  and  air  drag  on  a 
point  mass  non-lifting  reentry  body.  This  model  is  essentially  that 
of  Catalano;  the  derivations  given  here  are  somewhat  more  detailed, 
and  the  results  are  adapted  to  the  needs  of  programming  and  statistical 
investigations. 

In  the  following  derivation  vectors,  which  are  always  underlined,  are 
understood  to  be  expressed  in  column  matrix  form  in  a  particular  ortho¬ 
gonal  coordinate  system.  Thus,  the  physical  vector  from  the  radar  to 
the  tracked  object,  the  vector,  is  ^  in  the  unprimed  coordinate  system, 
in  the  primed  system,  and  JR"  in  the  double  primed  system.  These  are 
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all  diffeyent  matrices,  although  representing  the  same  pliysical  vector, 
and  have  a  common  length  S,  the  range  measurement.  When  the  physical 
vector  is  discussed,  the  particular  representation  H  will  be  used  for 
convenience.  These  conventions  will  be  followed  in  all  vector  notation, 
except  that  the  physical  unit  vectors  of  the  unprimed  system  (t^,  u^,  u^), 
of  the  primed  system  (Uj,  U2',  iij )»  of  the  double  primed  system 
—2.'‘  —3"^  represent  truly  different  vectors,  not  different  repre¬ 
sentations  of  the  same  vectors.  They  bear  the  prime  notation  to  distingmsh 
the  coordinate  set. 

Consider  then  an  oblate  earth  of  equatorial  radius  a^  and  eccentricity 
e,  rotating  at  the  siderial  rate  m.  Let  (»ij,  u^,  u.,)be  an  inertial  orthogonal 
coordinate  system  where  u^,  U2  are  in  the  equatorial  plane  and  u^  points 
toward  the  north  pole.  Asstune  a  radar  station  at  distance  a  from  earth's 
center,  at  geocentric  colatitude  <p  and  lying  at  longitude  0  with  respect  to 
the  u^  axis.  Figure  1  depicts  these  relations. 

"3 


Since  the  station  is  fixed  with  respect  to  the  earth, we  have 
9  =  u)  or  9  =  Gq  +  mt 
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The  polar  radius  of  the  earth  is  av/T-  . 

e 

Nov/  consider  a  cross-section  of  the  earth  in  the  plane  of  the  station. 
Figure  E.  h  is  the  altitude  of  the  station,  p  its  geodetic  latitude,  and  p* 
its  geocentric  latitude,  u^  is  the  equatorial  coordinate  in  the  station  plane. 


Figure  2. 

Also  shown  in  Figure  2  are  the  two  construction  circles  of  radii  a^  and 
a^\/ 1  -  e^  for  the  parametric  representation  of  the  elliptical  figure  with 
the  aid  of  the  angle  9 .  Thus,  the  ground  projection  of  the  station  has 
coordinates 

(2.1)  u  =  a  cos  P 
s  e 

/  2 

u,  =  a  V  1  -  e  sin  p  . 
i  e 
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a  cos  Ml 

(2.2)  Ug  =  -g  .  ^ 

yi  -  sin- 


+  h  cos  g, 

8 


(1  -  e^)  sin  g, 

u,  =  — -  +  h  sin  g, 


— 2 — —  s 

8in*g 


For  convenience  defire 


(2. 3)  aj 

^2 


/I  -  e-  sin*  g 
a  (l-e^) 

/i  -  e2  ain2  g  ® 


then  using  (2,  2),  the  station  vector  (coordinates)  in  the  (Uj,  U2,  u^)  system 
can  be  written 


(2.4)  a  = 


^2 


cos  g  cos  9 
cos  g  sin  9 
sin  g 


Note  that  a,  the  magnitude  of  ^  is  given  by 

CA  2  2  2  ^2.2 

(Z.  5)  a  =  aj^  cos  g  +  a2  sin  g 

« 

Noting  that  9  =  (u,  we  find  the  rates  of  change  of  ji  in  the  inertial  system 
to  be 


(2.6) 


U)  cos 


V- 


-  sin  0 
cos  6 
0 


cos 


cos  9 

sin  9 

0 
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The  station  observations  are  made  from  position  a  but  in  an  earth 
fixed  system  with  origin  at  a.  Define  an  orthogonal  system  at  ^  and 
make  observations  as  shown  in  Figure  3. 
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Figure  3. 


The  observed  object  is  at  point  O  at  position  vector  r  ;  the  observation 
vector  is  ^  .  In  the  local  system  (iij*,  u^’,  u^'),  u^'  is  perpendicular  to 
the  earth's  surface  (at  the  stations  ground  point),  points  tangentially 
east,  and  tangentially  north.  R  may  also  be  described  in  terms  of 
range,  azimuth,  and  elevation  relative  to  the  u'  system  as  shown  in 
Figure  4, 
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A  third  set  of  orthogonal  axes  may  then  be  constructed  with  u," 
along  with  the  R  direction,  U3"  in  the  A  direction  parallel  to  the  ground 
plane,  and  03"  in  the  E  direction  to  form  a  right-handed  system.  By 

projecting  the  u"  unit  vectors  on  to  the  u'  axes,  we  find  the  vector 
relations. 


(2.7) 

Hi” 

=  cos  E  sin  A  u^' 

+  cos  E  cos  A  U3' 

+  sin  E  U3' 

H2” 

=  -  sin  E  sin  A  u^' 

-  sin  E  cos  A  u^' 

+  cos  E  u^' 

-3" 

=  cos  A  Uj 

sin  A  u^* 

These  provide  an  orthogonal  transformation.  T^,  for  transforming  vector 

components  from  the  u"  to  the  u'  system. 


1 

cos  E  sin  A 

-  sin  E  sin  A 

cos  A  j 

cos  E  cos  A 

-  sin  E  cos  A 

-  sin  A  I  K 
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For  example,  the  physical  vector  ^  has  the  components  (R,  0,  0)  in 
the  u"  system.  In  the  u'  system,  the  coordinates  are 


(2. 9) 


R'  =  T.. 


R 

0 

0 


R  cos  E  sin  A 
R  cos  £  cos  A 
R  sin  E 


But  the  ^  vector  components  are  also  related  to  R '  by  an  orthogc-nal  matrix 
transformation. 


(2.10)  R  =  TR' 


From  the  geometry  implied  by  Figures  1  and  3,  we  can  ixxunediately  project 
unit  vectors  U2',  u^*  inertial  system  and  represent  them  as 

linear  combinations  of  unit  vectors  Uj,  U2,  system.  We  obtain 

the  vector  relations. 


-  sin  9  Uj^  +  cos  0U2 

-  sin  jj,  cos  0  Uj  -  sin  sin  ®  £2  1^3 

cos  \i  cos  0]^  +  cos  p,  sin  6  U2  ®^®  —s 


These  relations  immediately  yield  the  required  transformation. 


(2.11)  T 


sin  9  -  sin  p  cos  0 

cos  0  -  sin  p  sin  0 
0  cos  p 


cos  p  cos  9 
cos  p  sin  9 
sin  p 


3.  Equations  of  Motion 

From  Figure  3,  we  see  that 


(2.12)  r  =  a  +  R 
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Then  using  (<J.  13)  and  differ entiating,  we  find  the  succession  of  relations 

(2.13)  r  =  a  +  T  R' 

r  =  k  +  TR'  +  TR' 

Y  =  k*  +  T  R '  +  21  R '  +  T  R ' 

By  Newton's  law  (r  is  in  an  inertial  system)  r  may  be  equated  to  the  specific 
force  vector  applied  to  the  reentry  body. 


(2.14)  k*  +  T  R=  +  2T  R'  +  T  R'  =  F/m 

Now  by  direct  differentions,  the  left  side  of  (2.14)  is  evaluated  using  (2.6), 
(2.  S)i  (2.  9),  and  (2.11).  First  the  T  derivatives  may  be  evaluated  in  the 
following  form. 


(2.15)  T  = 


=  Tu) 


u)  cos  9 
0)  sin  9 
0 

0 

sin  p 
-  cos  p 


u)  sin  p  sin  9 
>  (u  sin  p  cos  6 
0 

-  sin  p  cos  p 

0  0 

0  0 


(u  cos  p  sin  9 
uu  cos  p  cos  9 
0 


D 


D 


(2.16)  T  =  T«») 


Toj‘ 


0 

sin  p 
-  cos  p 


1 

0 

0 


-  sin  p  cos  p 
0  0 

0  0 


.  2 
sin  p 

sin  p  cos  p 
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-  sin  p  cos  p 
2 

cos  p 


I 

f 

I 
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In  order  to  evaluate  R'  derivatives,  we  first  calculate  in  the  form 


I 

I 


1 

I 

I 

i 

I 


(2.17)  T.  = 


Then 


-EsinEsinA+AcosEcosA  -EcosEsinA-AsinEcosA  -AsinA 
•  •  •  •  * 
-EsinEcosAoAcosEsinA  -EcosEcosA+AsinEsinA  -AcosA 

•  • 

EcosE  -EsinE  0 


=  T. 


0  -  E  -  A  cos  E 

•  • 

E  0  A  sin  E 

•  • 

A  cos  E  .-A  sin  E  0 


=  T 


and 

(2.19)  R' 


=  X. 


=  X. 


R 

RE 

RA  cos  E 

•  • 

R 

RE  +  RE  1  +  T... 

•  •  •  •  •  • 

RA  cosE  +  RA  cos  E  -  RAE  sin  E 


••  *  7  *  7  y 

R  -  RE  -  RA  cos  E 
RE  +  2RE  +  RA^sin  E  cos  E 


RA  cos  E  -t  2RA  cos  E  -  3RAE  sin  E 


R 

RE 

RA  cos  E 


u 


By  directly  substituting  and  rearranging  these  results,  the  left  side 
of  (2.14)  takes  the  form 


I 


e  — “ 

R 

r 

R 

1 

(2.18)  R'  = 

0 

+ 

0 

1 

0 

0 

1 
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4,  Forces 

We  turn  now  to  a  consideration  of  the  epecific  force  on  the  right  side 
of  (2.14).  We  shall  include  two  forces,  gravity  and  air  drag.  The  gravita¬ 
tional  specific  force  upon  a  body  at  O  in  Figure  s  can  be  expressed  in  terms 
of  an  infinite  series.  Explicity  writing  terms  only  out  to  the  first  zonal 
harmonic  (J2)»  series  becomes 


(2.23)  F  /m  = 


r  1  +  I 
—  1  ^ 


(r/a^)‘ 


(1-5  cos  tf...  )  +  . . .  1  + 


coscp*  3  — ^ 
r  (r/a^) 


+  .  . . 


Here  £  is  the  body  position  vector  from  earth's  center,  r  is  its  magnitude, 
and  is  the  angle  it  sub'sndr  with  the  polar-axis.  Gj^  and  are  gravita¬ 
tional  constants  and  u^  is  a  unit  vector  in  the  north  polar  direction.  In 
using  (2.  23)  further,  we  shall  retain  only  the  terms  shown.  Applying 
(2.  4),  (2.  9),  (2, 13)  and  writing  out  the  vector  components  in  (2.  23)  in 
the  inertial  system,  we  obtain 


(2.  24)  F  /m  = 

O 


-■f  (1-5cosSj 


+  TT..  0 


cos  9^^  3 


(r/a^)' 


a^  cosM'  cos  6 
a^  COSM'  sin 6 


a2  sin^ 


Finally,  to  bring  out  a  uniform  TT^  factor  on  the  left  side,  we  note  that 
TT^  is  a  unit  matrix  and  terms  lacking  this  factor  need  to  be  multiplied 

bv  tI  r^.  Thus  (2.  24)  can  be  written  in  the  form 
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(2.25)  F  /m  =  -^  1+4 

-*  7^  ^  (r/a^)‘ 


(1-5  cos  cp^ ) 


(a2  -  «ij^)  sinpicosp,  + 

2  X  -2 

cos  pi  +  a^sin  p. 


cos  cp^  3 


(r/a^)' 


TT.,  <  T:  cos  ji 


sin  p.  . 


Now  introduce  the  abbreviations 


2  2 

(2.  26)  c^  =  cos  ^  +  a^  sin  p. 


c^  =  (aj^  -  a2)  sinp.  cosp, 


and  perform  the  T.I^  multiplications  to  obtain 


(2.27)  F  /m  = 

O 


r  3  j- 
L  (r/i 


(1-5  cos‘ 


%  ) 


R  +  Cj^  sin  E  -  C2  cos  E  cos  A 
TXj.  Cj^  cos  E  *-2  *1“  ^  ^  + 

C2  sin  A 


cos  «p. 


(r/a^) 


cosp  cos£cosA+  sinpsinE 
TTjj,  sinpcosE  -  coe  ;•  ?inEcosA 
-  cospsinA 


I 
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The  use  of  this  result  requires  an  evaluation  of  r  and  cos  We 

note  from  (2. 12)  and  (2. 13)  that  r  =  a  +  R  =a  +  TR'  so  that 


(2.28)  r  = 


t 

r  r 


=  a 


^  +  R^  +  2a^  T  R' 


=  +  R^  +  2R  [0,  -  Cj] 


cos  £  sin  A 
cos  £  cos  A 
sin  £ 


2  2 

=  a  +  R  +  2R  (Cj^  sin  £  -  c^  cos  £  cos  A) 


Also 


(2.  29)  cos  cp^ 


t 


=  i  (^+TR')‘u3 

=  —  I.  a2  Sinn  +  R  (cos  cos  £co8A  +  sin  t^sinE)] 

The  drag  specific  force  on  a  body  at  altitude  D  above  the  earth, 
traveling  with  velocity  V  relative  to  the  earth  is  given 

(2.30)  F^m  =  -|g  p(D)  «  VV 

where  gP  (D)  is  an  empirical  function,  the  air  density  in  weight  units  at 
altitude  D,  and  »  is  a  drag  coefficient  which  characterizes  the  body. 
Since  R'  is  the  relative  velocity  vector  expressed  in  the  earth  fixed  u* 
system,  we  can  write  V  in  the  inertial  system  as 

• 

R 


(2.31)  V  =  T  R'  =  TT^  I  RE 

PA  cos  £ 


b 
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an  approximation  to  D.  This  approximation  is  very  good  for  small  e. 
As  with  the  similar  expressions  (2.1),  the  point  O,.  can  be  located  in 
the  parametric  form 

(2.  34)  =  a^  cos 

u,  =  a  >/  1  -  e^  sin 

j  e  , 

But  also 

(2.35)  Uj.  =  (r-d)  sin  <P.., 

=  (r-d)  cos  9... 


Eliminating  P^.,  from  (2.  34),  (2,  35),  we  find 

2 


=  (r-d)' 


(si 


2 

sin  cp 


cos  cp 


and,  therefore, 
(2.  36)  D  ^  d 


%n/  ^ 

/ — rri! — 

v  1  -  e  sin  cp 

•v 


5.  State  Model  Equations 

Expressions  (2.  27)  and  (2.  33),  when  added,  complete  the  right  side  of 
the  differential  equations  of  motion  (2.14).  We  have  already  evaluated  the 
left  side  in  (2.  22).  Note  thac  all  these  expressions  have  a  common  matrix 
factor  TTo.  which  we  can  discard  or  eliminate  by  multiplying  throughout  by 
This  final  step  in  effect  transforms  all  vectors  into  the  u”  system 
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which  is  an  orthogonalized  radar  coordinate  (R,  E,  A)  system.  Thus,  the 
component  differential  dynamical  equations  for  the  motion  of  the  reentry 
body  in  radar  polar  coordinates  can  be  read  off. 

In  order  to  conform  to  requirements  of  the  Kalman  model,  that  the 
dynamical  equations  be  first  order  differential  equations  as  in  (1. 1),  we 
identify  as  dynamical  variables  both  R,  E,  A  and  R,  E,  A.  Since  we  shall 
be  interested  in  estimating  the  drag  parameter  «  as  well,  it  will  also  be 
included  in  the  state  vector  x  .  Thus,  define 


(2.  37) 


X 


R 

A 

E 

• 

R 

A 

E 

ot 


where  the  orde.*  of  variables  E,  A  has  been  interchanged  to  match  a  current 
convention.  Then  the  first  three  component  eqviations  of  state  (1. 1)  are 
actually  kinematical. 


(2.  38) 


dR 

■HF 


dt 


d£ 

■ar 


while  the  second  three  are  obtained  from  the  derived  equations  of  motion  as 
noted  above: 
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(2.  39) 


(2. 40) 


Range  equation 


dR  _  „  •  2  .  „  12  2  „ 

-  RE  +  RA  COB  E 


-  2u)R[£  cos^  8inA  + AcosE  (siniicoBE  -  cosiJ>sinEco8A)3  + 

+  u)  R  [l  -  (cospicosEcoBA  +  siniisin  E)^]  + 

2 

+  ajU)  coBp.  (coBjiBinE  -  BinptcosEcosA)  + 


M 


\ 


(1-5  COB  9 


j.)  [R  +  Cj  si 


sinE-c  cosEcobA]  + 
2 


cobEcobA  +  sin  sinE  ]  + 


^gp(D)a  VR 


Asiznuth  eqviation  (divide  by  RcosE) 


dt 


2RA 


+  2AE  tan  E 


~  cosE  I  R^  lisinEcosA-  sin  (icoaE)  +  E  (sin  ^l8inE+co8lJ/  co8Eco8A)j 

2 


uu 


+  COB  n«inA  (BinjisinE  +  cos  (icosEcosA) 

,  2 
a.  u) 

■*’RH5iE  »inp,coB»i8inA 


M 


r  RcosE 

2G, 


M 


cosE 


(1-5co8^9jJ 


c^sinA 


cos^sinA 


-  ^  gP  (D)  «  V  A 
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(2.41)  Elevation  equation  (Divide  by  R) 


dE  2RE  12  .  „  „ 

" - —  -  A  sinEcosE  + 


dt 


+  2uJ 


R 


cospiainA-i-AcosE  (siniisinE+cosixcosEcosA) 


2 


+  (u  (8inp.8inE+co8  (icosEcoeA)  (cob  ixsinEcosA-sinixcosE)  + 


*1  2 

+  (u  co8^(8iniJi8inEco8A  +  cosI^cobE)  -f 


[i+l  hCr)  a-5co.\)j 


2G 


M 


COB 


r  R 


<P« 


5  COB  (f^)| 

2. 


c^  cosE+c^  sinEcosA 


a . 


j 


■ [“ 


sin^cosE-coBp.  BinEcoBAI  + 


■] 


-  ^  g  p(D)  Of  V  E 


Ab  a  final  component  equation  of  the  Byatem  model  (1. 1),  we  append  the 
trivial  relation 


(2.42)  ^  =  0 


ThuB  for  the  present  we  aeaume  the  drag  parameter  to  be  constant. 


6.  Observation  Model  Equations 

The  next  step  in  constructing  a  model  for  recursive  trajectory 

estimation  is  to  delineate  the  form  of  the  observations  (1.  2),  Since  our 

•  •  • 

dynamical  description  is  in  an  R.  E,  A,  R,  E,  A  system,  which  is  directly 
the  set  of  observed  quantities,  the  relations  (1.  2)  have  a  very  simple  form. 
To  be  specific,  let  vts  suppose  that  at  each  observation  time  we  measure 
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R,  A,  E,  R.  Then  the  obeervation  equation  components  are 


(2.43) 

II 

R  +  Wj 

*2  " 

A  +  Wj 

*3  = 

E  +  Wj 

*4  " 

R  +  w . 

4 

where 

Wi,  Wy 

Wy  and  are  the  respective  additive 

measurements 

with  a  known  covariance  matrix  W. 

We  have  noted  previously  that  statistical  calculations,  the  actual 
process  of  track  smoothing  and  prediction,  require  a  linearization  of  both 
the  dynamical  equations  (2.  38)  through  (2,  42)  and  the  observation  equations 
(2. 43).  But  the  latter  are  already  in  linear  form;  the  H  matrix  (1. 12)  for 
this  set  of  measurements  then  has  the  simple  form, 

1  0  0  0  0  0  0 

0  1  0  0  0  0  0 

0  0  1  0  0  0  0 

0  0  0  1  0  0  0 


(2.44) 

H  = 
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m.  SPHERICAL  EARTH  CASE  -  JACOBIAN  MATRIX 
1.  SiTnpIifiod  Equations  of  Motion 

While  the  observation  matrix  H  is  trivial,  the  J  matrix  (1.8)  for  our 
dynamical  system  is  very  complicated.  We  have  described  in  section  13 
the  use  of  the  J  matrix  in  the  tracking  procedure.  It  is  used  in  (1.  14)  to 
extrapolate  the  E  matrix  between  data  times  while  the  dynamical  equations 
(1.  13)  are  used  to  extrapolate  the  estimated  trajectory.  Essentially,  then, 
the  J  matrix  is  only  used  in  error  propagation,  and  hence  need  not  be  cal¬ 
culated  numerically  with  the  same  accuracy  as  the  f_  functions.  Tn  fact, 
the  only  effect  of  slightly  inaccurate  J's  is  to  produce  slightly  inaccurate 
E's  w'hich  result  in  less  than  optimal  estimates,  x's.  The  effect  is  anal¬ 
ogous  to  the  result  of  perturbing  the  value  of  a  variable  y  near  a  minimum 
value  of  a  function  g(y).  Relatively  large  excursions  of  y  are  often  tolerable 
before  the  value  of  g(y)  is  appreciably  raised. 

It  is  very  convenient  for  our  problem  to  simplify  J.  For  this  reason, 
and  with  the  above  justification,  we  first  reduce  our  f_  functions  (2.39), 
(2.40),  (2.41)  to  the  case  of  a  spherical  earth,  (e  =  0),  rxnd  drop  the 
gravitational  terms,  (J^  =  0),  These  simplifications  also  imply  that 


(3.1) 


ai  =  a2  a  =  Cl 


c^  =  0 


r^  =  a^  +  R^  +  2aRsinE 


coscpjj,  ~  sinp  +  R(cospcoe£cosA  +  sin|isinE)J 


D  =  r  -  a  +  h 

s 


Making  the  above  simplifications  we  note  that  (2.  38)  remains  unchanged 
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(3.2) 


dR 

dt 

dA 

dt 

dE 

dt 


R 


'  A 


=  E 


while  (2.39),  (2.40),  (2.41),  respectively  simplify  to: 


(3.3) 


(3.4) 


dR 

dT 


*  2  *2  2 
RE^  +  RA^cos'^E  + 


dA 

dt 


r  ■  "  1 

-  2u)R  j  EcosiisinA  +  AcosE(sinjiCosE  -  cos^iSinEcosA)  j  + 

2  r  21 

+  CD  -  (cos^iCOsEco8A  +  sinp,sinE)  j  + 

2 

+  au)  cosii(cosusinE  -  siniiCosEcosA)  + 

Gw  I 

-  -y  (R  +  asinE)  -  i  gp(D)erVR  . 


2RA 


+  ZAEtanE  + 


2(1)  I'R 


cosE 


\  R 

I  y  (cosp,sinEcosA  -  sinp.cosE)  + 

'  1 

+  E(sin(j,smE  +  cospcosEcosA)  j  + 


+  cosp.sinA(sin'^sinE  +  cosiicosEcosA)  + 

2 


-igP(D)QtVA  . 
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(-’•5)  - - — - A  sinEcosE  + 


1  R  ’  1 

+  2u)|^^  cosiisinA  +  Aco8E(8iniJisinE  +  cos^lCosEco8A) J  + 


4  'X  (sinjiSinE  4-  co8p,co8Eco8A)(co8nsinEcosA  -  sin^.co8E)  4- 
cl  3 

4-  ^  u)  cosii{sinu.sinEco8A  4-  co8p.co8E)  4- 
1 

— X —  acosE  -  ygp(D)aVE  . 
r  R 


and  again 


(3.6) 


der 

dt 


=  0 


2.  Jacobian  Matrix 


The  J  matrix  has  the  following  form  (whether  eimplified  or  not). 


(3.7) 


J  = 


S(fR»  fA»  ^E’ ^R’  ^A'^E' 


S(R.A, 

E,R. 

A.E, 

a) 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

0 

RR 

^kA 

■^RE 

T*  • 

‘'rr 

J*  ■ 

‘'ra 

T*  • 

‘'re 

•^Ra 

AR 

■^AA 

^AE 

T  •  • 

‘^AR 

T  **  • 

aa 

T  •  • 

‘'ae 

ER 

^EA 

^EE 

J  •  • 

‘'er 

J  •  • 

T  •  • 

‘'ee 

^EQf 

0 

0 

0 

0 

0 

0 

Performing  the  indicated  differentiona  with  the  right  aide  of  2)  through 
(3.6)  we  can  derive  the  individual  J  components.  We  note  the  expression 
for  V  in  (2.  32);  we  also  use  the  fact  that  for  any  state  variable,  x. 
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The  individual  simplified,  spherical  earth  J  matrix  components  are 
(3.9)  +  A^cob^F-  + 

-  Zuj^EcospsinA  +  Acos  E(sinpco8£  -  cospsinEcosA)^  + 

+  -  (cospcosEcosA  +  sinpsinE)^^  + 

-  — j-  +  — ^ —  (R  +  asinE)  + 

,  .  'z  :z  2„ 

-igp(D)0fRR  + 

"  jgp'(D)aR(R  +  asinE) • 


(3.10) 


(3.11) 


=  -  2u}Rcosp(EcosA  +  AsinEcosEsinA)  + 


+  2(1}  Rcospco8EsinA{cospcosEcosA  +  sinpsinE)  ‘t 


+  au)  cospsinpcosEsinA 


J^j,  =  -2RAsin£co8E  + 


+  2u)RA|^28inu,sinEcosE  +  co8p(co8^E  -  8in^E)co8Aj  + 

2 

-  2(1}  R(co8pco8£co8A  -(■  8inp8in£)(sinpco8E  -  co8p8inEco8A)  4- 


+  au}  co8p(co8pco8£  4-  8inp8lnEco8A)  4 
G^acosE 

,  aRco8E{R+a8inE)  + 
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(3.18) 


(3.19) 


(3.20) 


(3.21) 


(3.22) 


AWOOM 


£fe. 


1%  77.\  T  •  _  2RE  2ujR  _  .  ^ 

w*  ”  ”  2  “  — 2 —  cosjisinA  + 

R  R 


au) 


i. 

1 


m  h 


(3.24) 


Y~  cosn(sinp,sinEcosA  +  cosjicosE)  + 
M® 


R 

Gj^acosE  3G^aco8E(R  +  asinE) 


3  _2 
r  R 


r^R 


-  |^gp(D)tfER  + 


-  ^go*{D)aE{R  +  asinE)  — 

^  T 


R  *  2 

=  2u)C08-  cosA  -  Acos  EsinA)  + 


-  U) 


'cos4sinAr2cos|i,s  inEcosEcosA  -  8inu(cos^E  -  sin^E)]  + 


- U)^cosnsinp.8inEsuiA  , 


(3.26) 

(3.27) 


(3.28) 

(3.29) 


+  jgo(D)erE 


R^A^sinEcooE 


|•gp'(D)aEaRcosE-~ 


'ER 


^  cosjisinA  -  igp(D)erE 


=  -  2AsinEcosE  +  2u)CosE(sin|ASinE  +  cosiacosEcosA)  + 


•  R^Acos^E 


±gp(D)aE 


o  •  «. 


*  2  2 
E  R 


'EE 


-igp(D)V.igo(D)or^ 


=  -  J8P{D)EV 
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IV.  STATISTICAL  TESTING  PROCEDURES 


1.  Basic  Procedure 

It  is  evident  from  the  foregoing  description  that  the  equations  for 
tracking  in  radar  polar  coordinates  are  relatively  complex.  Although 
one  can  organize  the  calculations  to  take  advantage  of  common  trignono- 
metric  expressions,  etc. ,  in  the  formulas,  the  program  running  time 
for  such  calculations  considerably  exceeds  that  for  a  simpler  coordinate 
s/fatem  (earth  fixed  rectangular)  with  a  simpler  model  (spherical  earth). 

In  order  to  compare  the  radar  polar  coordinate  tracking  procedure  with 
simpler  versions  of  itself  and  with  other  procedures  a  systematic  testing 
method  is  desirable. 

Such  a  testing  method  can  be  implemented  by  a  computer  simulation 
which  statistically  compares  the  tracking  runs  made  on  a  given  trajectory 
(or  set  of  trajectories)  using  the  full  model  against  similar  runs  made 
with  a  simplified  or  alternate  model.  Let  us  call  the  full  model  the 
general  model  and  the  various  alternate  models,  special  models. 

The  statistical  comparisons  can  be  implemented  by  the  following  pro¬ 
cedure. 

a)  Select  a  given  trajectory  by  specifying  its  initial  conditions  and  the 
drag  characteristics  of  the  body.  Select  the  radar  site  location. 
Specify  statistical  a  priori  error  variances  for  the  initial  trajectory 
and  drag  (state)  variables.  Select  the  quantities  to  be  measured,  the 
data  sampling  rate,  and  the  measurement  error  covariance  matrix. 

b)  Pass  1.  Simulate  the  trajectory  of  the  reentry  object  by  integrating 
the  equations  of  motion  as  given  by  the  general  model  with  the  initial 
conditions  assumed  above.  This  generated  trajectory  is  assumed  to 
be  the  true  trajectory.  Record  its  values  x,j;(n)  (state  vector)  at  the 
sampling  times. 
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At  the  same  time  make  and  record  z(n).  These  are  measurements 
of  the  true  trajectory  with  an  added  random  noise  vector  chosen  by  a 
Monte  Carlo  process  from  a  multidimensional  normal  distribution  with 
the  assumed  covariance  matrix.  (Since  the  coordinate  system  of  current 
concern  is  radar-centered  polar,  the  H(n)  matrices  (2.44)  are  trivial  and 
need  not  be  specifically  computed  or  stored.  In  general,  however, 

Hjj.(n)  would  be  evaluated  on  the  reference  trajectory  and  recorded. ) 

Perform  recursive  tracking  as  outlined  in  section  I  using  again  the  full 
model  with  the  data  £(n)  and  the  above  input  parameters.  Evaluate  partial 
derivatives  (the  J  matrix  computation)  at  the  true  trajectory  points  rather 
than  at  the  estimated  trajectory  values.  This  tracking  procedure  provides 
an  ideal,  optimum  set  of  estimates  Xjjj(n|n-1),  x^j^  (njn)  and  their  associated 
covariances  I^(n|n)  which  are  recorded. 

The  recorded  quantities  provide  a  standard  against  which  alternate, 
simplified  tracking  procedures  may  be  compared. 

c)  Pass  2.  Using  pass  1  observations,  z(n),  as  data,  perform  the  track¬ 
ing  according  to  the  special  alternate  model  to  be  tested.  This  tracking 
generates  the  estimates  x(n|n-l),  x(n|n)  and  their  covariances  E(n|n-1), 
I(njn).  Here  we  do  not  assume  that  the  true  trajectory  is  known;  the 
J  matrix  is  computed  along  the  estimated  trajectory,  also  the  initial 
track  estimates  and  covariances  need  not  be  identical  to  those  of  pass 
1.  This  is  a  simulation  of  tracking  as  it  would  actually  be  performed 
in  the  real  world. 

During  the  pass  2  tracking  we  compute  a  series  of  comparison  statis¬ 
tics  to  be  described  below.  These  comparison  statistics  determine 
whether  the  special  model  tracking  has  been  significantly  degraded  from 
the  ideal. 

2.  Tests  of  Trajectory  Fit 

In  order  to  derive  these  statistics  it  is  convenient  to  view  the  Kalman 
tracking  process  as  a  very  general  method  for  least  squares  curve  fitting. 
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For  example,  suppose  that  a  curve  of  the  form 


(4.1) 


•  -an  ,  „  -bn  ,  r-  "Cn 
Ae  +  Be  +  Ce 


•with  parameters  A,  a,  B,  b,  C,  c  is  to  be  fitted  to  data  points,  z(n), 
n  =  1,  2,  3, ...  N.  This  problem  can  be  linearized  and  formulated  according 
to  the  Kalman  model  with  the  state  vector 


A 

a 

B 

b 

C 

c 


X  = 


dx 


having  the  trivial  dynamics  =  0.  The  computational  scheme  of  section 
I  can  be  immediately  adapted  to  this  problem. 

Suppose  we  are  sure  (the  null  hypothesis)  that  the  data  arise  from  a 
model  of  the  above  form;  then  we  may  ask  the  question  whether  the  ex¬ 
pression  (4.  1)  for  the  curve  can  be  simplified  by  dropping  the  third  term  , 
viz. , 


(4.2) 


.  -an  .  „  -bn 
Ae  +  Be 


A  rather  obvious  procedure  for  determining  the  suitability  of  the  special 
expression  (4.  2)  in  representing  the  data  is  to  perform  a  least  squares  fit 
of  (4.2)  to  the  data,  z(n),  and  to  examine  the  residuals 


(4.3) 


"(n)  -  (Ae”^°  -f  Be”^"  ) 


after  the  best  fit  has  been  obtained.  Under  the  null  hypothesis  the  expected 
size  of  the  residuals  can  be  determined.  If  the  residuals  are  large,  in  some 
collective  sense,  with  respect  to  the  e}q>ected  magnitude  of  the  discrepancy. 
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then  the  special  model  {4.  2)  is  detectably  poorer  than  the  general  model 
(4.  1)  in  its  ability  to  fit  the  data.  The  simplification  has  significantly 
degraded  the  model. 

This  inspection  of  the  residuals  can  be  used  in  a  similar  fashion  to 
test  alternatives  to  (4.  1)  even  though  they  may  have  a  completely  differ¬ 
ent  form.  For  instance,  instead  of  (4.  1)  we  might  try  a  polynomial  repre¬ 
sentation 

(4.4)  An^  +  Bn  +  C 

to  fit  the  data.  Under  some  circumstances  we  might  find  that  residuals 
with  (4.4)  are  satisfactorily  small  and  that,  therefore,  a  polynomial  curve 
will  also  fit  the  data.  Note  that  such  a  type  of  residual  test  does  not  deter- 

^  A  A  ^  A  A  A 

mine  whether  or  not  the  best  fit  parameters  A,  B,  C  of  (4.  4)  or  A,  a,  B,  b, 
or  (4.  2)  are  close  to  or  correspond  in  any  sense  to  the  best  fit  parameters 

A  A  A  A 

A,  a,  B,  b,  C,  c  of  the  general  model  (4.  1).  Indeed,  while  the  general 
model  may  have  a  physical  basis  for  its  form,  the  model  (4.  4)  may  be 
entirely  empirical.  Yet  (4.4)  may  equally  well  fit  the  data;  its  estimates 

A  ^  A 

A,  B,  C,  then,  are  all  that  are  needed  to  represent  a  fitted  curve.  In 
effect,  they  extract  all  of  the  available  information  from  the  curve  data. 

The  remainder,  the  residuals,  appear  typically  random. 

The  above  test  of  fit,  then,  is  basically  an  answer  to  the  question  of 
whether  or  not  simplifications  or  alterations  in  the  curve  fitting  or,  in 
general,  the  tracking  model  produce  a  loss  of  information  during  the  pro¬ 
cess  of  estimation.  A  negative  result  indicates  that  the  reduced  model 
complexity  is  unaccompanied  by  the  penalty  of  information  loss. 

In  order  to  formalize  the  above  ideas  let  us  assume  that  our  model 
has  been  linearized  and  that  its  statistics  are  gaussian  (measurement  noise, 
random  driving  terms  in  state  equation,  initial  conditions).  Then,  given 
the  model,  one  can  write  down  the  joint  density  distribution  of  the  data 
^(1),  £(2), .  . , ,  z(N)  given  the  initial  state  value  Xq.  This  distribution  has 
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the  joint  normal  form 

(4.5)  p{z(l),£(2), . . . ,  z{N)*Xq)  =  Ce 


1q 


where  C  is  a  constant  and  Q  is  a  quadratic  form  in  the  z's.  By  the  chain 
rule  of  conditioning  we  can  also  write  (4.  5)  in  the  form 

(4.6)  p(£{l),  z(2), . .  .  ,£(N)j  Xq)  + 

=  p(£{l/I  Xq)p{£(2)1  £(1),  Zq)p(z(3)I  z{l),  z(2),Xq)  + 

...  p(z(N)!z(l),....z(N-i),x  )  + 


Ce 


-  i(Qj  +  Q2  >  ...  +  Qj^) 


where  is  the  quadratic  form  in  the  density  distribution  p(z{n)\^l), 

.  .  .  ,  z{n-l),z  ). 

A 

It  is  now  simple  to  identify  £(nln-l).  the  expected  value  of  z_(n)  given 
the  previous  data  and  initial  conditions.  Since  ^(n)  =  H(n)x(n)  +  w(n), 

(4.7)  z(nin-l)  =  H(n)x(nln-1)  , 

for  the  noise  w{n),  being  independent  of  previous  data,  has  conditional 
expectation  zero.  Thus  has  the  form 

(4.8)  =  |^£(n)  -  H(n)x{nln-i)^  S"^|'£(n)  -  H(n)x(n!n-l)J 

where  iw  the  conditional  covariance  matrix  of  £(r). 

By  direct  calculation  we  find 
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(4.9)  =  cov|^£{n)  -  {'i(n)x(nln- 

=  cov|^H(n){}c(n)  -  x(nln-l))  +  ^(n)^ 

=  H(n)£{n!n-l)H^(n)  +  W(n)  , 

since  x{n)  -  x(n|n-l)  and  w^(n)  are  independent.  Thus  reviewing  (4.5),  (4.6), 
(4.8),  (4.9)  we  see  that  the  data  distribution  (4.  5)  which  may  also  be  termed 
the  likelihood  function  uses  the  quadratic  form 

N  .  -1 

(4.  10)  Q  =  I  (z(n)  -  H(n)x(nin-1))'  H(n)I(nln-l)H'^(n)  +  W(n)  1  • 

n=l  ^ 

•  (£(n)  -  H(r)x{nln-1))  . 

This  quadratic  form,  being  in  the  exponent  of  a  joint  normal  distribution, 
has  a  chi-square  distribution.  It  is  evident  that  the  vector  difference 

A 

£(n)  -  H(n)x(n|n-1)  is  a  type  of  residual  showing  the  deviation  of  the  i-ew  data 
£(n)  from  its  position  expected  on  the  basis  of  the  past  data  and  initial  con¬ 
ditions.  Thus  expression  (4.  10)  is  a  measure  of  data  fit  over  the  span 
n  =  1,  Z, .  - .  ,N.  Large  values  of  Q  indicate  bad  fit,  small  values  indicate 
good  fit.  Statisticians  will  also  note  that  Q  is  -2  times  the  natural  log  of 
the  likelihood  function  and  that  the  value  of  this  ftmetion  indicates  the  degree 
to  which  the  data  conforms  to  the  assumed  mod’l. 

The  dimensionality  of  Q,  the  number  of  degrees  of  freedom,  is,  on  the 
face  of  it,  Np,  where  each  z  vector  has  p  components.  This  will  be  true  in 
nondegenerate  cases  where  the  residuals  in  (4.  1C)  and  their  associated 
covariances  are  not  infinitely  small  or  mfinitely  large  and  are  consist  vtly 
computed.  The  one  important  practical  variation  to  this  rule  is  the  fre¬ 
quently  useful,  but  inconsistent,  assumption  of  an  initial  vector  having 
finite  (small)  components  but  initial  T(l|0)  matrix  with  iminite  (very  large) 
variance  components.  This  inconsistency  manifests  itself  in  the  computation 
of  Q  where  some  of  the  initial  residuals  are  much  smaller  than  their  purported 
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variance.  Thus,  no  preceptible  contribution  to  Q  results.  In  order  to 
correct  for  this  effect,  we  must  subtract  from  Np  the  number  of 
vector  components  so  handled.  The  result  is  the  correct  number  of 
degrees  of  freedom,  u ,  for  Q.  For  example,  in  a  tracking  model  where 
four  observations  are  taken  at  each  time  step  and  Xq  has  seven  compo¬ 
nents,  the  two  extreme  instances  are  as  follows. 


(^.11) 


r 


4N 


jj4N  -  7 


initial  x,  Z  consistent 

initial  x  components  all  «  than  initial  Z  variances 


3.  Variations  on  the  Test  of  Fit 

The  Q  statistic  may  now  be  used  in  the  two  pass  testing  procedure 
mentioned  previously.  £(n)  are  generated  on  the  1st  pass  with  the  full 
model;  x{nln-l)  on  the  2nd  pass  with  the  special  model  to  be  tested.  We 
have  two  choices  for  the  covariance  matrices  Zj^^Cnln-l)  from  pass  1  or  the 
estimated  covariances  on  pass  2.  The  pass  1  covariances  actually  supply 
the  required  statistic  because  the  pass  2  estimated  ccvariances  are  subject 
to  additional  computational  errors  that  the  estimated  trajectory  must  be 
substituted  for  the  true  trajectory  in  the  linearization  process.  (In  our 
current  application  H{n)  is  not  affected  by  this  distinction  since  in  radar 
polar  coordinates  no  linearization  of  the  observational  equations  is  re¬ 
quired.  In  general,  however,  H(u)  would  be  derived  by  a  linearization  and 
a  corresponding  choice  of  evaluation  on  the  true  trajectory,  —  from  pass 
1  —  or  on  the  estimated  trajectory  would  be  afforded. 

For  certain  special  purposes  we  might  prefer,  nevertheless,  to  use 

the  estimated  covariances  to  calculate  Q.  Our  interest  is  dvte  to  the  fact 

t 

that  the  matrix  H£H  +  W  and  the  residual  vector  z  -  Hx  occur  in  the 
estimate  and  covariance  recxirsions  (1.  15),  (1  16)  of  the  tracking  process 
as  well  as  in  the  Q  terms.  Numerical  difficulties  can  occur  in  the  tracking 
which  lead  to  insta-bilities  if  the  estimated  matrix  HZH^  +  W  departs  signi¬ 
ficantly  from  the  ideal.  If  such  difficulties  are  present,  the  residual  vector 
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will  not  compare  well  with  H  IH*’  f  W  and  the  Q  term  will  tend  to  be  large. 
Note  that  here  the  situation  is  essintially  the  reverse  of  the  previous  one 
where  H  +  W  is  known  to  be  correct  and  the  residuals  are  compared 
with  it.  In  the  present  situation  the  accuracy  of  the  estimated  HEK^  +  W 
is  also  in  question.  For  such  a  situation  the  Q  statistic  provides  a  useful 
measure  of  the  discrepancy  although  no  specific  numerical  probabilities 
can  be  attached  to  the  measure  because  the  Q  sampling  distribution  is  then 
not  simply  "ihi-square.  By  comparing  the  Q's  obtained  by  using  both  ideal 
and  estimated  E‘s  we  can  qualitatively  assess  the  additional  covariance 
estimate  errors  and  stability  properties  of  the  tracking  recursions. 

We  note  from  expression  (4.  10)  that  a  Q  value  may  be  computed  for 
any  particular  number  of  successive  data  points  N.  Strictly  speaking,  a 
numerically  correct  test  of  fit  requires  us  to  fix  N,  (e.g. ,  to  include  ail 
the  data)  before  calculating  Q  since  the  successive  increasing  values  of  Q 
are  not  statistically  independent.  It  is  convenient,  nevertheless,  to  monitor 
tliese  successive  values  so  as  to  get  an  approximate  idea  of  the  way  the 
trajectory  fit  is  progressing.  In  order  to  expedite  this  monitoring  it  is  con¬ 
venient  to  convert  Q  to  an  equivalent  unit  normally  distributed  variable,  X. 
Thus  we  map  the  chi- squared  distribution  and  regions  shown  in  Figure  6a 
onto  the  unit  normal  distribution  and  regions  in  Figure  6b. 


f/Q) 
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Thjs  mapping  depends  on  tht  number  of  degrees  of  freedom  and  is  given 

(5) 

approximately  by  the  formula 

H.!2)  X  .  it  -  . 

VT/(9yl 

On  the  basis  of  chance  we  expect  that  X  will  vary  roughly  between 
-2  and  +2.  If  X  >  2,  then  the  residuals  are  too  large,  the  fit  is  bad  and 
the  model  is  suspect.  If  X  <  -2,  the  fit  is  too  good;  this  circumstance  is 
usually  due  to  a  programming  erroi , 

4.  Test  of  Estimate  Accuracy 

Now  suppose  that  the  special  model  fit  has  been  tested  and  is  satis¬ 
factory.  A  further  question  is  whether  the  estimates,  x{nln),  produced 
by  the  special  model  are  close  to  the  true  trajectory  values  £Cjj.(n),  That 
this  need  not  be  ths  case  is  shown  by  our  curve  fitting  example  noted  pre¬ 
viously.  The  state  variables  describing  the  polynomial  (4.4)  of  the  special 
model  are  not  directly  comparaVle  to  the  variables  in  the  exponential 
expression  (4.  1)  dewcribing  the  general  model;  yet  the  fit  can  be  excellent. 

In  many  instancec,  however,  the  intent  of  the  model  specialization  will 
be  not  only  to  preserve  the  fit  but  also  to  preserve  the  significance  of  the 
model  parameters.  In  these  instances  we  hope  to  find  close  correspon¬ 
dence  between  x^(n)  components  and  all  or  some  of  the  x{n!n)  components. 

An  appropriate  test  for  this  correspondence  (for  the  entire  vector)  is 
to  form  the  estimate  error,  x(n|n)  -  x,),(n),  and  compare  it  with  its  theoreti¬ 
cal  ideal  covariance  ^^^.(nln)  in  a  statistic, 

(4.  13)  P  =  j^x(nln)  -  x,^(n)J  L"^(n|n)  |^x(nln)  -  £,^(n)J  . 

This  statistic  can  be  computed  for  each  n  and,  except  for  the  firat  few 
sample  points  where  degeneracies  may  occur  as  noted  earlier,  it  has  a 
chi-.<jquare  distribution  with  q  degrees  of  freedom,  where  q  is  the  dimension 
of  x^. 
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If  x(n'n)  is  significantly  distorted  fronn  the  true  the  P  statistic 

will  be  large,  A  typical  situation  for  a  two  dimensional  state  vector  as 
shown  in  Figure  7. 


Figure  7. 

V/e  see  an  error  ellipse  of  constant  P  and  constant  error  probability  as 
determined  by  rjjj(n|n).  Let  us  assume  that  this  ellipse  marks  the  critical 
P  level.  If  the  error  vector  lies  outside  this  ellipse  the  error  is  signifi¬ 
cantly  large.  Such  an  error  vector  e  ,  is  shown  in  the  figure. 

However,  note  vector  e^j^.  This  vector  also  is  significantly  large. 

But  owing  to  the  high  statistical  correlation  between  error  components, 
causing  a  very  skewed  error  ellipse,  the  individual  components  of 
when  tested  one  at  a  time  against  their  variances,  are  not  individually 
large.  Thus  it  is  possible  for  the  combined  error  vector  test  to  show  up 
significant  distortions  even  though  individual  components  of  error  are  well 
within  their  respective  expected  limits. 

For  this  reason  it  is  appropriate  to  examine  specific  components  or 
subsets  of  the  vector  x(n|n/  for  the  desired  correspondences  based  on  the 
physical  requirements  of  the  problem.  For  example,  in  a  trajectory  track¬ 
ing  problem  we  might  separate  the  trajectory  position  and  velocity  components 
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of  X  into  OK^:  subvector  and  the  drag  parameter  components  into  another 
and  test  each  separately  against  corresponding  components  of  x^j^. 

One  further  test  variation  is  useful.  Under  some  circumstances  our 
special  model  parameters  may  be  chosen  to  be  specific  functions  of  the 
general  model  parameters.  Cr  we  may  be  interested  in  a  subsidiary  esti¬ 
mation  of  suck  parameter  functions  although  the  recursive  estimation  does 
not  produce  them  directly.  For  example,  let  us  consider  a  vector  y  sf 
functions  of  the  general  model  state  at  time  n. 

(4.  14)  _jr{n)  =  ^{x^(n))  . 

Suppose  that  our  estimation  process  produces  the  estimates  v(nln)  by  the 
calculation, 

(4.  15)  X(n|n)  =  _^(x(nin))  . 

* 

Then  the  errors  ^(n)  -  y;{nln)  may  be  compared  with  the  theoretical  co- 
variance  as  before  in  a  chi-square  test.  For  this  conxparison  we  need 

(4.  16)  cov|^^{nln)^  =  cov|^jy{x{n|n)J 

r  6  V  "1  ...  r 

“  [-Sr]  «v<*(nln))}_— J 


or 

(4.  17)  S{n)  =  r{K)  E^(n!n)  F^n) 

where  r(n)  is  the  partial  derivative  matrix  of  the  v  functions  with  respect 
to  the  X  variables  evaluated  along  the  true  trajectory.  If  we  assume  no 
degeneracy,  then  the  niimber  of  degrees  of  freedom  for  the  chi-squared 
test  statistic. 
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(4.  17)  P  =  j^2(n)  -  i(nln)]  S"^n)r_jr(n)  -  X'"!")] 
is  equal  to  the  nuiTiber  of  x  components. 
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NOTATION  AND  SYMBOL  LIST 


General  Notation 


matrix  transpose 


matrix  inverse 


vcvtor 


vector  components 


matrix  components 


probability  density  distribution 


dt  '  ' 


time  derivative,  equation  of  state 


time  derivative,  dynamic  equations  of  motion 


quantity  evaluated  on  reference  trajectory 


deviation  from  the  reference  value 


X  (or  other  quantity)  at  time  i 


cov(  ) 


covariance  operator 


estimate 


x(nlk)  1 

^nlk)  J 


estimate  of  x(n)  and  covariance  of  estimate 
given  data  up  to  and  including  time  tj^ 


section  IV  notation  for  covariances,  estimates 
and  observations  calculated  on  pass  1 


A.  a,  B,  b,  C,  c 


general  parameters  in  section  IV  discussion 
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I 

I 

i 
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Coordinate  Systems  (orthogonal) 

Uj,  U2,  Uj  inertial 

unit  vectors  Uj,  2i3  (Figure  1) 

Uj^',  U2',  Uj  radar  centered,  topographic 

unit  vectors  Uj,  U2',  U3'  (Figure  3) 

xij",  U2",  Uj"  radar  centered,  along  line  of  sight 

unit  vectors  1^",  ^2"*  —3  *  (Figure  4) 

Conventions 


Let  Tf  be  a  physical  vector.  Then  we  express 


Y 

its  u^,  U2>  coordinates 

Xl 

its  u^',  U2',  u^  coordinates 

Y  *  * 

its  Uj",  U2",  uy  coordinates 

Y  = 

111 

=  iri  =  11"  1 

p  m 


Transforixiatioii  s 


Y'  =  Y" 


Y  =  T  Y' 


S,nnibol8 


a,  a 


a,,  a 


V  “2 


radar  station  vector,  a  =  |a| 
station  parameters 
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Symbols  (cont.  ) 


equatorial  earth  radius 


radar  azimuth  coordinate 


drag  parameter 


B  B 


c, ,  c 


1'  ^2 


auxiliary  angles  in  ellipse  geometry 


derived  quantities  dependent  on  station  parameters 


constant  in  probability  distribution 


altitude  of  tracked  object  above  earth, 
approximation  to  D 


eccentricity  of  earth 


®1'  ®2»  ~a'  ®b 


estimate  error  components  and  vectors 


radar  elevatron  coordinate 


functions  on  right  side  of  equations  of  state 


F,  F  ,  F , 
— ’  — g'  — d 


vector  force,  gravity,  drag 


principal  gravitational  constant 


arbitrary  functions  of  state  variables 


-  Y* 

partial  derivative  matrix  of  ^  ,  F. .  =  $ — - 

M  a*. 

observation  functions 
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Symbols  (cont. ) 
h 

s 

J 

■^2 

k 

m 

IJI>  u' 

n 

N 

y 


o»  O  j.i 


p 


q 

Q 

£ 

R.  B 


station  height 

jacobian  matrix,  J..  = -r - 

•'  '  IJ  a  Xj 

first  zonal  harmonic  gravitational  constant 

auxiliary  constant 

mass  of  tracked  object 

geodetic,  geocentric  latitude 

index  denoting  t^ 

total  number  of  observations 

number  of  degrees  of  freedom  for  chi-square 

position  points 

dimension  of  jt 

geocentric  colatitude;  station,  object 

dimension  of  x 

quadratic  form,  test  of  fit 

object  vector  from  earth  center 

object  vector  from  radar,  radar  range  coordinate 
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Symbols  (cont. ) 
gp  (D) 

S 

t 

T.X. 

•Y* 

e 

u 

y,  y 

w 

W 

•r 

W.. 

X 

X 

1 

z 

tt) 


air  density  (weigh?  units)  at  altitude  D 
covariance  of  z 

covariance  x  estimates 

time 

transformation  (rotation)  matrices 
station  longitude 

(with  subscripts  and  primes)  various  coordinate 
system  unit  vectors  and  continents 

velocity,  speed  of  object  relative  to  earth 

observation  noise 

observation  noise  covariance 
state  equation  driving  term 
driving  term  covariance 
state  vector 

normal  variate  equivalent  to  chi-square  variate 
variables  equal  to  functions  of  x 
observation  vector 
earth  rotation  rate  (siderial) 
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